package com.alex.statistics.method.spatialAnalysis;

import com.alex.statistics.pojo.request.spatialAnalysis.GlobalMoranIRequest;
import com.alex.statistics.pojo.result.spatialAnalysis.GlobalMoranIResult;
import org.geotools.feature.FeatureCollection;
import org.geotools.geojson.feature.FeatureJSON;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.geom.Point;

import org.opengis.feature.simple.SimpleFeature;
import org.springframework.stereotype.Service;


import org.opengis.feature.Property;

import java.io.StringReader;
import java.util.*;
import java.util.stream.Collectors;

@Service
public class GlobalMoranIService {

    private final FeatureJSON featureJson = new FeatureJSON();

    /**
     * 通过GeoJSON计算全局Moran's I空间自相关指数
     */
    public GlobalMoranIResult calculate(GlobalMoranIRequest request) {
        long startTime = System.currentTimeMillis();
        validateGeoJsonInput(request);

        // 1. 解析GeoJSON获取空间要素
        List<SimpleFeature> features = parseGeoJsonFeatures(request.getGeoJson());

        // 2. 提取属性值（value）
        List<Double> values = extractValuesFromFeatures(features, request.getValueField());
        int n = values.size();

        // 3. 生成空间权重矩阵
        double[][] weights = generateSpatialWeights(features, request);

        // 4. 计算基础统计量
        double mean = values.stream().mapToDouble(Double::doubleValue).average().orElse(0);
        double sumSq = values.stream().mapToDouble(v -> Math.pow(v - mean, 2)).sum();
        double s0 = Arrays.stream(weights).flatMapToDouble(Arrays::stream).sum();

        // 5. 计算全局Moran's I指数
        double moranI = computeMoranI(values, weights, n, mean, sumSq, s0);

        // 6. 构建结果
        GlobalMoranIResult result = new GlobalMoranIResult();
        result.setMoranIndex(moranI);

        // 7. 计算额外统计量（如果需要）
        if (request.getIncludeStats()) {
            double expected = -1.0 / (n - 1);
            double variance = computeVariance(n, s0, sumSq, weights);
            result.setExpectedValue(expected);
            result.setVariance(variance);
            result.setZScore((moranI - expected) / Math.sqrt(variance));
        }

        // 8. 蒙特卡洛显著性检验
        if (request.getPermutations() > 0) {
            double pValue = computePValue(values, weights, moranI, request.getPermutations());
            result.setPValue(pValue);
            result.setSignificant(pValue < 0.05);
        }

        // 9. 判断空间模式
        result.setPattern(interpretPattern(moranI, result.getPValue()));
        result.setComputeTime(System.currentTimeMillis() - startTime);

        return result;
    }

    /**
     * 解析GeoJSON为要素列表（支持FeatureCollection）
     */
    private List<SimpleFeature> parseGeoJsonFeatures(String geoJson) {
        try {
            FeatureCollection featureCollection= featureJson.readFeatureCollection(new StringReader(geoJson));

            List<SimpleFeature> features = new ArrayList<>();
            // 获取迭代器并遍历
            try (var iterator = featureCollection.features()) {
                while (iterator.hasNext()) {
                    features.add((SimpleFeature)iterator.next());
                }
            }

            return features;
        } catch (Exception e) {
            throw new IllegalArgumentException("GeoJSON解析失败：" + e.getMessage(), e);
        }
    }

    /**
     * 从要素中提取指定字段的数值
     */
    private List<Double> extractValuesFromFeatures(List<SimpleFeature> features, String valueField) {
        return features.stream()
                .map(feature -> {
                    Property property = feature.getProperty(valueField);
                    if (property == null) {
                        throw new IllegalArgumentException("GeoJSON中未找到属性字段：" + valueField);
                    }
                    Object value = property.getValue();
                    if (value instanceof Number) {
                        return ((Number) value).doubleValue();
                    } else {
                        throw new IllegalArgumentException("属性字段'" + valueField + "'不是数值类型");
                    }
                })
                .collect(Collectors.toList());
    }

    /**
     * 根据请求参数生成空间权重矩阵
     */
    private double[][] generateSpatialWeights(List<SimpleFeature> features, GlobalMoranIRequest request) {
        int n = features.size();
        List<Point> centroids = features.stream()
                .map(feature -> ((Geometry) feature.getDefaultGeometry()).getCentroid())
                .collect(Collectors.toList());

        switch (request.getWeightType()) {
            case KNN:
                return generateKnnWeights(centroids, request.getKnnK());
            case DISTANCE_THRESHOLD:
                return generateDistanceThresholdWeights(centroids, request.getDistanceThreshold());
            case QUEEN:
                return generateQueenWeights(features);
            case ROOK:
                return generateRookWeights(features);
            default:
                throw new IllegalArgumentException("不支持的权重类型：" + request.getWeightType());
        }
    }

    /**
     * K最近邻权重矩阵（行标准化）
     */
    private double[][] generateKnnWeights(List<Point> centroids, int k) {
        int n = centroids.size();
        double[][] weights = new double[n][n];

        for (int i = 0; i < n; i++) {
            Point pointI = centroids.get(i);
            // 计算与其他点的距离并排序
            Map<Integer, Double> distanceMap = new HashMap<>();
            for (int j = 0; j < n; j++) {
                if (i != j) {
                    distanceMap.put(j, pointI.distance(centroids.get(j)));
                }
            }
            // 取最近的K个点
            List<Integer> nearestNeighbors = distanceMap.entrySet().stream()
                    .sorted(Map.Entry.comparingByValue())
                    .limit(k)
                    .map(Map.Entry::getKey)
                    .collect(Collectors.toList());
            // 行标准化（每个邻居权重为1/K）
            double weight = 1.0 / k;
            for (int j : nearestNeighbors) {
                weights[i][j] = weight;
            }
        }
        return weights;
    }

    /**
     * 距离阈值权重矩阵（行标准化）
     */
    private double[][] generateDistanceThresholdWeights(List<Point> centroids, double threshold) {
        int n = centroids.size();
        double[][] weights = new double[n][n];

        for (int i = 0; i < n; i++) {
            Point pointI = centroids.get(i);
            List<Integer> neighbors = new ArrayList<>();
            // 收集距离≤阈值的点
            for (int j = 0; j < n; j++) {
                if (i != j && pointI.distance(centroids.get(j)) <= threshold) {
                    neighbors.add(j);
                }
            }
            // 行标准化
            if (!neighbors.isEmpty()) {
                double weight = 1.0 / neighbors.size();
                for (int j : neighbors) {
                    weights[i][j] = weight;
                }
            }
        }
        return weights;
    }

    /**
     * Queen邻接权重矩阵（面状要素，行标准化）
     */
    private double[][] generateQueenWeights(List<SimpleFeature> features) {
        int n = features.size();
        double[][] weights = new double[n][n];

        for (int i = 0; i < n; i++) {
            Geometry geomI = (Geometry) features.get(i).getDefaultGeometry();
            List<Integer> neighbors = new ArrayList<>();
            // Queen邻接：共享边界或顶点
            for (int j = 0; j < n; j++) {
                if (i != j) {
                    Geometry geomJ = (Geometry) features.get(j).getDefaultGeometry();
                    if (geomI.touches(geomJ) || geomI.intersection(geomJ).getDimension() >= 0) {
                        neighbors.add(j);
                    }
                }
            }
            // 行标准化
            if (!neighbors.isEmpty()) {
                double weight = 1.0 / neighbors.size();
                for (int j : neighbors) {
                    weights[i][j] = weight;
                }
            }
        }
        return weights;
    }

    /**
     * Rook邻接权重矩阵（面状要素，行标准化）
     */
    private double[][] generateRookWeights(List<SimpleFeature> features) {
        int n = features.size();
        double[][] weights = new double[n][n];

        for (int i = 0; i < n; i++) {
            Geometry geomI = (Geometry) features.get(i).getDefaultGeometry();
            List<Integer> neighbors = new ArrayList<>();
            // Rook邻接：仅共享边界（维度为1的交集）
            for (int j = 0; j < n; j++) {
                if (i != j) {
                    Geometry geomJ = (Geometry) features.get(j).getDefaultGeometry();
                    Geometry intersection = geomI.intersection(geomJ);
                    if (intersection.getDimension() == 1) { // 交集为线（边界）
                        neighbors.add(j);
                    }
                }
            }
            // 行标准化
            if (!neighbors.isEmpty()) {
                double weight = 1.0 / neighbors.size();
                for (int j : neighbors) {
                    weights[i][j] = weight;
                }
            }
        }
        return weights;
    }

    /**
     * 计算Moran's I指数
     */
    private double computeMoranI(List<Double> values, double[][] weights, int n, double mean, double sumSq, double s0) {
        double numerator = 0;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                if (i != j) {
                    numerator += weights[i][j] * (values.get(i) - mean) * (values.get(j) - mean);
                }
            }
        }
        return (n / s0) * (numerator / sumSq);
    }

    /**
     * 计算理论方差
     */
    private double computeVariance(int n, double s0, double sumSq, double[][] weights) {
        double s1 = 0.5 * Arrays.stream(weights).flatMapToDouble(Arrays::stream).map(w -> Math.pow(w + w, 2)).sum();
        double s2 = Arrays.stream(weights).mapToDouble(row -> Math.pow(Arrays.stream(row).sum(), 2)).sum();

        double nSq = (double) n * n;
        double term1 = (nSq - 3 * n + 3) * s1 - n * s2 + 3 * s0 * s0;
        double term2 = (nSq - n) * s1 - 2 * n * s2 + 6 * s0 * s0;
        double term3 = (nSq - 3) * s0 * s0;

        return (term1 / ((n - 1) * (n - 2) * (n - 3) * s0 * s0))
                - (term2 / ((n - 1) * (n - 1) * (n - 2) * n * s0 * s0))
                + (term3 / ((n - 1) * (n - 1) * n * s0 * s0));
    }

    /**
     * 蒙特卡洛置换检验计算P值
     */
    private double computePValue(List<Double> values, double[][] weights, double observedI, int permutations) {
        int countExtreme = 0;
        Random random = new Random();
        int n = values.size();
        double mean = values.stream().mapToDouble(Double::doubleValue).average().orElse(0);
        double sumSq = values.stream().mapToDouble(v -> Math.pow(v - mean, 2)).sum();
        double s0 = Arrays.stream(weights).flatMapToDouble(Arrays::stream).sum();

        for (int p = 0; p < permutations; p++) {
            // 随机置换值
            List<Double> shuffled = new ArrayList<>(values);
            Collections.shuffle(shuffled, random);

            // 计算置换后的I值
            double randomI = computeMoranI(shuffled, weights, n, mean, sumSq, s0);
            if (Math.abs(randomI) >= Math.abs(observedI)) {
                countExtreme++;
            }
        }

        return (countExtreme + 1.0) / (permutations + 1.0); // 避免P=0
    }

    /**
     * 解释空间模式
     */
    private String interpretPattern(double moranI, Double pValue) {
        if (pValue == null || pValue >= 0.05) {
            return "空间分布随机（无显著自相关）";
        }
        return moranI > 0 ? "空间正相关（高值与高值聚集，低值与低值聚集）" : "空间负相关（高值与低值相邻，分布离散）";
    }

    /**
     * 验证GeoJSON输入参数
     */
    private void validateGeoJsonInput(GlobalMoranIRequest request) {
        if (request.getGeoJson() == null || request.getGeoJson().trim().isEmpty()) {
            throw new IllegalArgumentException("GeoJSON不能为空");
        }
        if (request.getValueField() == null || request.getValueField().trim().isEmpty()) {
            throw new IllegalArgumentException("必须指定valueField（属性字段名）");
        }
        if (request.getWeightType() == null) {
            throw new IllegalArgumentException("必须指定weightType（权重矩阵类型）");
        }

        // 权重类型对应的参数验证
        if (request.getWeightType() == GlobalMoranIRequest.WeightType.KNN) {
            if (request.getKnnK() == null || request.getKnnK() <= 0) {
                throw new IllegalArgumentException("KNN权重类型必须指定正整数的knnK");
            }
        } else if (request.getWeightType() == GlobalMoranIRequest.WeightType.DISTANCE_THRESHOLD) {
            if (request.getDistanceThreshold() == null || request.getDistanceThreshold() <= 0) {
                throw new IllegalArgumentException("距离阈值权重类型必须指定正数的distanceThreshold");
            }
        }

        // 置换次数验证
        if (request.getPermutations() < 0) {
            throw new IllegalArgumentException("置换次数不能为负数");
        }
    }
}
